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ABSTRACT 

We investigate the time evolution and spatial variation of the stellar initial 
mass function (IMF) in star-forming disk galaxies by using chemo dynamical sim- 
ulations with an IMF model depending both on local densities and metallicities 
([Fe/H]) of the interstellar medium (ISM). We find that the slope (a) of a power- 
law IMF {N{m) oc m~") for stellar masses larger than IMq evolves from the 
canonical Salpeter IMF {a ~ 2.35) to be moderately top-heavy one {a ~ 1.9) 
in the simulated disk galaxies with starbursts triggered by galaxy interaction. 
We also find that a in star-forming regions correlates with star formation rate 
densities (Esfr in units of Mq yr~^ kpc~^). Feedback effects of Type la and II 
supernovae are found to prevent IMFs from being too top-heavy {a < 1.5). The 
simulation predicts a ~ 0.23 log Ssfr + 1-7 for logSsFR > —2 (i.e., more top- 
heavy in higher Ssfr); which is reasonably consistent well with corresponding 
recent observational results. The present study also predicts that inner regions 
of starburst disk galaxies have smaller a thus are more top-heavy {da/dR ~ 0.07 
kpc~^ for i? < 5 kpc). The predicted radial a gradient can be tested against 
future observational studies of the a variation in star-forming galaxies. 

Subject headings: galaxies: ISM — galaxies: stellar content — stars: formation 
— stars: luminosity function, mass function 



1. Introduction 

A growing number of recent observational studies have suggested that stellar initial 
mass functions (IMFs) can be different in galaxies with different star-formation activities, 
dynamical parameters (e.g., stellar velocity dispersion), environments, and redshifts (e.g., 
Hoversten & Glazebrook 2008; van Dokkum 2008; Meurer et al. 2009, M09; Treu et al. 
2010; Cappellari et al. 2012; Ferreras et al. 2012; van Dokkum & Conroy 2012). One 
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of important questions relating to possibly variable IMFs is as to whether IMFs can be 
top-heavy in starburst galaxies (e.g., Elmegreen 2009 for a recent review). It has been still 
controversial whether the IMF is almost universal in different galaxies or depends on local 
environments such as densities and metallicities of star-forming gas clouds (e.g., Bastian et 
al. 2010; Kroupa et al. 2012). 

Recent observational studies on the physical relations between the slopes of power-law 
IMFs {a, where a = 2.35 is the canonical Salpeter IMF) and galaxy properties (e.g., mass 
densities and SFRs) have provided new clues to the nature of IMFs in actively star-forming 
galaxies. M09 investigated the extinction-corrected flux ratio of Hq and far- ultraviolet (FUV) 
light and found that the flux ratios are different between different galaxies and depend on 
surface densities of Hq, flux (Sh„) and it!-band luminosity densities of galaxies (Sr). Their 
observations are consistent with a decreasing (i.e., becoming more top-heavy) with increasing 
Eh„ and Er. Gunawardhana et al. (2011, Gil) investigated a of ~ 33000 galaxies and their 
correlations with SFRs (in units of Mq yr~^) and SFR densities (in units of Mq yr~^ kpc~^) 
and found correlations between a and SFRs (a ~ —0.36 log(SFR) -|- 2.6) and between a and 
SsFR (ct ~ — 0.3 log(SsFR) + 1.7). It should be noted here that they take a linear scaling 
between Ha luminosity and the SFR based on a = 2.35 in deriving the SFR. In reality, a 
linear scaling does not apply if a is varying: the exact scaling dependent on the behavior 
of the IMF in detail (e.g. the mass limits and any changes in a at the low mass end) is 
beyond the scope of this paper to determine. Although these observed correlations appear 
to strongly suggest top-heavy IMFs in starburst galaxies, the origin of the correlations has 
been discussed by only a few theoretical models in a quantitative manner (Pflamm-Altenburg 
et al. 2009; Weidner et al. 2011; Kroupa et al. 2012). 

Pipino & Matteucci (2004) showed that the observed colors of massive elliptical galaxies 
can be reproduced better by their chemical and photometric evolution models with slightly 
flatter (i.e., top-heavy) IMFs. Nagashima et al. (2005) demonstrated that the observed 
chemical abundances (e.g., [Mg/Fe]) of elliptical galaxies can be explained by top-heavy 
IMFs. These previous works, however, did not investigate how IMFs evolve in the starburst 
phases of elliptical galaxy formation. Although starbursts have been demonstrated to be 
triggered by galaxy interaction in previous numerical simulations (e.g., Noguchi & Ishibashi 
1986), the evolution of IMFs during starbursts has not been investigated in the simulations 
of galaxy interaction. Thus it is theoretically unclear whether or not the IMFs of galaxies 
become top-heavy when they are in the starburst phases. 

The two main purposes of this Letter are as follows. Firstly, we discuss whether or not 
the IMFs in starburst galaxies can become top-heavy by using our original chemodynamical 
simulations that incorporate the new physical IMF model proposed by Marks et al. (2012, 
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M12). Secondly, we test whether or not the present model can reproduce the observed 
correlations between a, SFRs, and Ssfr derived by Gil. We investigate a at individual 
star-forming regions in disk galaxies at each time step so that we can investigate not only 
the time evolution of IMF but also local variation of a. We focus our investigation on the 
IMF evolution of galaxies with starbursts triggered by strong galaxy-galaxy interaction and 
merging. In the present paper, we compare between the observed and simulated a — Ssfr 
relations, even though the observed relation would be less precisely derived by assuming the 
fixed Salpter IMF for the estimation of SFRs. 

2. The model 

In order to perform numerical simulations of galaxies on GPU clusters, we use our 
original chemodynamical code (Bekki 2013) that is a revised version of our previous code 
('GRAPE-SPH'; Bekki 2009). The new GPU- version chemodynamical code allows us to 
investigate the time evolution of chemical abundances and dust properties of galaxies self- 
consistently. We adopt a disk galaxy model that broadly mimics the structural and kine- 
matical properties of baryonic and dark matter components of the Milky Way Galaxy (as 
listed in Binncy & Trcmaine 2007). The disk galaxy is composed of dark matter halo, bulge, 
stellar disk, and gas disk. The density distribution of the dark halo is represented by the 
'NFW' profile predicted by the cold dark matter cosmology (Navarro et al. 1996), and the 
total mass, the virial radius, and the c-parameter are set to be lO^^M©, 245 kpc, and 10, 
respectively. 

The stellar bulge is assumed to have the Hernquist profile with the total mass of 10^^ Mq 
and the disk size of 3.5 kpc. The radial (R) and vertical (Z) density profiles of the stellar 
disk are assumed to be proportional to exp(— 7?/7?o) with scale length Rq = 3.5 kpc and 
to sech^(Z/Zo) with scale height Zq — 0.7 kpc, respectively. In addition to the rotational 
velocity caused by the gravitational field of disk, bulge, and dark halo components, the initial 
radial and azimuthal velocity dispersions are assigned to the disc component according to 
the epicyclic theory with Toomre's parameter Q — 1.5. The maximum circular velocity of 
the disk is 222 km s~^ for the adopted structure parameters of the stellar and dark matter 
components. 

The gas disk represented by SPH particles has a gas mass fraction of 0.1 and an expo- 
nential radial profile with the scale-length of 3.5 kpc and the size of 35 kpc. The gas disk has 
an initial temperature of lO'' K and a metallicity gradient with the central [Fe/H] of 0.34 and 
the slope {d[Fe/}i]/dR) of —0.04 dex kpc^^ (Andrievsky et al. 2004). The gas particles are 
converted into new stars if the following physical conditions are met: (i) the local dynamical 
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time scale is shorter than the sound crossing time scale (mimicking the Jeans instability) 
and (ii) the local velocity field is identified as being consistent with gravitationally collapsing 
(i.e., div v< 0). The feedback energy of 10^^ erg is imparted to neighboring gas particles 
around a new star when the star becomes a Type la or Type II supernova (SN la and SN 
II, respectively). 

For a new star formed from Jeans-unstable gas, the IMF slope (a) for stellar masses 
(mg) larger than IMq is investigated. In the following, the high-mass IMF slope (originally 
0:3 in the Kroupa IMF; Kroupa 2001) is referred to as a for convenience. M12 proposed that 
a depends on mass densities (pci) and [Fe/H] of proto-cluster gas clouds as follows: 

a = 0.0572 X [Fe/H] - 0.4072 x log,,{j^^^^—^) + 1.9283 (1) 

This equation holds for Xtji > —0.87, where Xth = — 0.1405[Fe/H] -|- logio( io6jt^^'pc-a ); and 
a = 2.3 for xth < -0.87 (M12). We slightly modify the M12's IMF model in the following 
two points. Firstly, we do not introduce the threshold Xth, mainly because the a-dcpcndence 
at a; < Xth. is not so clear (not so flat as M12 showed) owing to a small number of data points. 

Secondly, we estimate pci from local gas density pg of Jeans-unstable gas particles in 
equation (1). This is because the present simulation can not resolve the rather high-density 
cores of star-forming or cluster-forming molecular gas clouds. The local value of Pci at each 
Jeans-unstable gas is estimated by multiplying pg by a constant kg. 

Pel = hPs- (2) 

The derivation process of kg is as follows. The Jeans-unstable gas is assumed to be star- 
forming giant molecular clouds and thus have the following size-mass scaling relation derived 
from the observed mass-density relation by Larson (1981): 

Since the equation (1) is based largely on the observed properties of the Galactic globular 
clusters (GCs hereafter) and nearby star clusters, a reasonable kg can be the typical density 
ratio of GCs to GC-host GMCs. We here consider that (i) pci should correspond to a typical 
mean mass density for GCs, (ii) typical GC mass (Mg^) and size (Rgc) arc 2 x IO^Mq and 
3pc, respectively (Binney & Tremaine 2007), (iii) original GCs just after their formation 
from GMCs should be ~ 10 times more massive than the present ones (e.g., Decressin et al. 
2010; Marks & Kroupa 2010; Bekki 2011), and (iv) a star formation efficiency of GC-host 
GMCs is ~ 0.1. For Mg^^ = 2 x IO^Mq and i^g^c = 283 pc in the typical GC-host GMC, a 
reasonable kg {— {Rgmc/ Rgc)^) is estimated to be 8.4 x 10^. The present models do not show 
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large a (~ 3; very steep IMF) that is needed to explain the observed ifa/FUV flux ratios 
of low surface brightness galaxies (LSBs) in M09. This is mainly because star formation 
can occur only in higher density gaseous regions where a can be smaller in the present 
models. Although the original IMF model by M12 is theoretically derived from observations 
on physical properties of GCs and star clusters, we assume that the IMF model applies for 
all new stars in each bin (i.e., not just for star clusters). 

The disk galaxy is assumed to interact with a companion galaxy represented by a point- 
mass particle. The mass-ratio of the companion to the disk galaxy is a free parameter denoted 
as 1712- The initial distance of the two galaxies is fixed at 280 kpc and the orbital eccentricity 
(e) and the orbital pericenter (rp) are free parameters. The angle between the z axis and the 
vector of the angular momentum of the disk is denoted as 9. We mainly describe the results 
of the flducial model with m.2 — 1, e — 1 (i.e., parabolic), — 35 kpc, and ^ = 0° (i.e., 
prograde interaction). In order to discuss briefly the parameter dependences, we investigate 
the following flve representative models: weaker (m2 = 0.3) and stronger {ni2 = 3) tidal 
interaction, retrograde interaction [9 — 180°), low surface brightness (LSB; Rq = 5.5 kpc), 
and no SN feedback effects. For comparison, we also investigate the isolated model (i.e., no 
tidal interaction) and the major merger one in which two identical disks merge with each 
other in a prograde-prograde manner. For the merger model, rp = 17.5 kpc and e = 0.8 are 
adopted. The total number of particles is 716700 for isolated and interaction models and 
1433400 for the merger model. The gravitational softening length is 2.3 kpc for dark matter 
and 0.25 kpc for baryonic components. In the following, T in a simulation represents the 
time that has elapsed since the simulation started. 



3. Results 

Figure 1 shows how a changes during a starburst triggered by galaxy interaction in 
the flducial model. After the pericenter passage of the companion galaxy, the strong tidal 
fleld can compress the gas disk so that rather high-density gaseous regions can be formed. 
As a result of this, the star formation rate (SFR) of the disk galaxy can be signiflcantly 

enhanced (up to ~ IOA/q) and the mean a (a) increases from ~ 2.2 (pre-starburst) to ~ 1.9 
(during starburst). The time evolution of SFR is clearly synchronized with that of a, which 
implies that a can be described as a function of SFR in galaxies. Individual local regions of 
the disk galaxy during the triggered starburst have different a ranging from 1.5 to 2.6 
owing to locally different gas densities and [Fe/H]. These results thus demonstrate that IMFs 
can become moderately top-heavy in starburst galaxies because of the formation of rather 
high-density and Jeans-unstable gaseous regions. 



-6- 



Figure 2 shows that a in individual radial bins (i.e., azimuthally averaged a) of the 
disk galaxy correlates with SFR densities (Ssfr) in the fiducial model. Although dispersions 
are large at a given Esfr, radial bins with higher SgpR are more likely to show smaller 
a (i.e., more top-heavy IMFs). These simulated positive correlations between a and Ssfr 
are qualitatively consistent with the corresponding observational results by Gil. Figure 
3 describes the radial a gradients of the star-forming disk at pre-starburst (T=0.3 Gyr), 
strong star burst (T=l.l Gyr), and post-star burst phase (T=1.7 Gyr). Clearly, more top- 
heavy IMFs (i.e., smaller a) can be seen in the inner regions of the interacting disk galaxy 
for the three epochs. The disk has da/dR of 0.07 kpc~^ for the central 5 kpc at T = 1.1 and 
1.7 Gyr. The derived positive a gradients can be seen in all models of the present study and 
thus regarded as a robust prediction. 

Figure 4 shows the a — Esfr correlation of individual radial bins with log Esfr > —2 
in the five different interacting galaxies with starbursts. Clearly, radial bins with higher 
Ssfr can have more top-heavy IMFs (i.e., smaller a). A least-squares fit to these simulated 
star-forming regions gives 

Q; = -0.231ogSsFR + 1.7, (4) 

which is similar to the observed relation of a ~ — 0.31ogEsFR + 1-7 in Gil. It should be 
noted here, however, that the observed relation is derived from averaged a and Ssfr (over 
individual galaxies) whereas the simulated points arc derived from the average of star forming 
regions in each radial bin, with the points belonging to different galaxies. It should be also 
noted that the simulated a — Ssfr relations can be slightly different between different models. 
For example, a — —0.28 log Ssfr + 1.7 in the fiducial model, which is pretty close to the 
observed relation by Gil, whereas a — —0.25 log Esfr+1-7 in the model with m2 = 0.3. The 
LSB model can show more top-heavy IMFs in some regions, because high-density gaseous 
regions can form owing to efficient gas-transfer to the central region during tidal interaction 
hence original LSB galaxies achieves a high surface brightness in the interaction. 

Figure 5 shows that (i) number distributions of a in the fiducial model and the model 
without SN feedback effects are bimodal (i.e., two peaks) and (ii) they are significantly 
different between the two models. The first peak at lower a ('low-a peak') is due to a 
starburst and the second peak at higher a ('high-a peak') represents the mean a for pre- 
starburst IMFs. The model without SN feedback has the low-a peak at lower a (a ~ 1.4), 
which means that the IMF is more top-heavy in comparison with the fiducial model. This 
difference between the two models clearly demonstrates that the feedback effects of SN la 
and SN II can suppress the formation of overly high-density, Jeans-unstable gaseous regions 
and therefore prevent IMFs from being too top-heavy {a < 1.5). Figure 5 also shows that 
the location of the low-a peak in the merger model is coincident with that of the fiducial one, 
which confirms that the IMF can become moderately top-heavy in major galaxy merging. 
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4. Discussion and conclusions 

The consistency between the observed and simulated a — Ssfr relations implies that 
the adopted IMF model (Ml 2) can be useful in simulating the time evolution and spatial 
variation of a in galaxies. However, this does not necessarily mean that the adopted IMF 
model for Jeans-unstable gas clouds is the best one that can be adopted for discussing 
galaxy evolution caused by varying IMFs. Larson (1998) adopted an IMF that is similar 
to the Salpeter IMF at the upper end and flatter one below a characteristic stellar mass 
(mch, which is typically ~ 1-/^© in the present day universe) and discussed a number of 
observational results in terms of time- varying mch- Following the IMF model proposed 
by Larson (1998), Dave (2008) and Narayanan & Dave (2012) considered IMFs with 
depending on redshifts and thereby discussed the origin of the observed redshift-evolution of 
the relation between stellar masses and SFRs in galaxies. It should be noted here that the 
IMF model by M12 is observationaUy supported whereas the mch variation model by Larson 
is a working hypothesis for which there is not much observational evidence. 

Given that a growing number of observations have recently revealed possible evidence 
for varying IMFs (e.g., Kroupa et al. 2012), it would be important for theoretical studies of 
galaxy formation and evolution to incorporate varying IMFs in their models. It is however 
unclear what varying IMF models should be adopted for galaxy formation and evolution 
studies. The results of numerical simulations on galaxy formation would depend on whether 
IMF models with varying mch (with fixed a) or those with varying a are adopted. A 
reasonable varying IMF model would need to explain recent observational results on IMFs in 
nearby galaxies such as correlations between /Jc^/FUV flux ratios and galaxy properties (e.g., 
M09) and between H^^ equivalent width and optical colors (e.g., Hoversten & Glazebrook 
2008; Gil) in a self-consistent manner. 

One robust prediction in the present study is that inner regions of actively star-forming 
disk galaxies have smaller a (i.e., more top-heavy). Such radial gradients of IMFs are 
also discussed in the context of the integrated galaxy IMF theory (Pflamm-Altenburg & 
Kroupa 2008). The predicted positive a gradient is a natural consequence of the adopted 
IMF model: the metallicity gradient act to flatten the IMF with increasing radius, but 
the pg-dependence is stronger so the new gradient is a steeper IMF at larger radius. Since 
observational studies have not yet extensively investigated radial (or azimuthal) gradients of 
IMFs in nearby galaxies, it is not clear whether the predicted a gradient is qualitatively or 
consistent with observations. We suggest that observed radial gradients of IMFs and their 
correlations with local galaxy properties (e.g., local gas densities and i?-band local surface 
luminosity densities) would give some constraints on whether IMF models with varying a or 
mch should be adopted for galaxy formation studies. 
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Recent observations have suggested that (i) massive early-type galaxies have the steeper 
(i.e., bottom-heavy) IMF in the mass range O.IMq to IMq (e.g., van Dokkum & Conroy 2011) 
and (ii) the IMF slope depends on stellar velocity dispersions of the galaxies (e.g., Ferreras et 
al. 2012). Since metal-rich star formation can lead to "bottom-heavy" IMFs for rris < IMq 
in M12, it is our future study to incorporate an IMF for rris < IMq that depends on physical 
properties of ISM in our chemodynamical simulations for elliptical galaxy formation. Such 
more sophisticated chemodynamical simulations with variable IMF slopes below and above 
rus — IMq will enable us to discuss the origin of the observed possibly bottom-heavy IMFs 
in massive early- type galaxies (e.g., Cenarro et al. 2003). 

We are grateful to the anonymous referee for constructive and useful comments. K.B. 
acknowledges the financial support of the Australian Research Council throughout the course 
of this work. Numerical computations reported here were carried out both on the GPU 
clusters (Pleiades and Fornax) at the University of Western Australia. 
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Fig. 1. — Time evolution of the distance of two interacting galaxies (top), the star formation 
rate (SFR; middle), and the IMF slope a (bottom) in the fiducial model. The big red and 
blue filled circles are the mean a for all new stars at each time interval in the fiducial 
(interaction) model and for those at the solar neighborhood (6 < i? < 10 kpc) in the isolated 
model, respectively. The solar neighborhood in the isolated model has almost constant a 
(~ 2.3 — 2.4; close to the Salpeter IMF), which strongly suggests that the adopted varying 
IMF model does a good job in predicting the IMF evolution. 
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Fig. 2. — The locations of star-forming regions on the logSgpR — a plane in the fiducial 
model. For convenience, —a is plotted against logEspR- Ssfr are estimated at 10 radial 
bins {R < 17.5 kpc) for star-forming regions (i.e., for only new stellar particles with ages 
less than ~ 10^ yr) at each time step of the simulation. 
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Fig. 3. — Radial profiles of —a for i? < 10 kpc at three epochs, pre-starburst (T = 0.3 Gyr; 
solid blue), during strong starburst (T = 1.1 Gyr; dotted red), and post-starburst (T = 1.7 
Gyr; dashed green) in the fiducial model. The dotted black line indicates the mean a for 
all new stars formed during the simulation. The radial gradients are estimated by using all 
new stellar particles with a. The central 5 kpc at T = 1.1 and 1.7 Gyr shows da/dR ~ 0.07 
kpc~^. 
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Fig. 4. — The locations of individual star-forming regions on the log Ssfr — a plane. Dif- 
ferent colors shows star-forming regions in different models: fiducial (blue), weaker inter- 
action with 1712 = 0.3 (red), stronger interaction with m2 = 3 (green), retrograde orbital 
configuration (magenta), and LSB model (cyan). The observed a — Ssfr relation by Gil 
(a ^ — 0.3 log(SsFR) + 1-7) and the best-fit simulated one for all of the star-forming regions 
are shown by dotted and dashed black lines, respectively. 
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Fig. 5. — The number distributions of a for the fiducial model (solid blue), the model 
with no SN feedback effects (dotted red), and the merger model (dashed green). The a 
distributions are derived by using a of all new stars in the models. 



